

A NEW LIE SERIES METHOD FOR THE 
NUMERICAL INTEGRATION OF ORDINARY 
DIFFERENTIAL EQUATIONS, WITH 
AN APPLICATION TO THE RESTRICTED 
PROBLEM OF THREE BODIES 


by Siegfried Filippi 

George C. Marshall Space Flight Center 
Huntsville, Ala. 






NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


WASHINGTON, D. C. 


APRIL 1967 



NASA TN D-3857 


A NEW LIE SERIES METHOD FOR THE NUMERICAL 
INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS, 
WITH AN APPLICATION TO THE RESTRICTED 
PROBLEM OF THREE BODIES 
By Siegfried Filippi 

George C. Marshall Space Flight Center 
Huntsville, Ala. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 - CFSTI price $3.00 



TABLE OF CONTENTS 


Page 

SUMMARY i 

I. INTRODUCTION 2 

n. GROBNER’S METHOD OF LIE SERIES 2 

A. General 2 

B. Proof of Convergence 4 

C. The Lie Series for the Solution of Ordinary 

Differential Equations 12 

D. Generalized Lie Series 16 

HI. APPLICATION TO THE RESTRICTED PROBLEM 

OF THREE BODIES 17 

A. General 17 

B. The Restricted Problem of Three Bodies . 18 

REFERENCES 34 


BIBLIOGRAPHY 


35 



A NEW LIE SERIES METHOD FOR THE NUMERICAL 
INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS 
WITH AN APPLICATION TO THE RESTRICTED 
PROBLEM OF THREE BODIES 1 


By 

Siegfried Filippi 2 


SUMMARY 


In this paper a new method for the numerical solution of problems of in- 
itial value in ordinary differential equations, originally developed by W. Grobner, 
will be described. The first step will be to describe the most important theorems 
of W. GrSbner’s Lie series with reference to the generalized Lie series of Filatov 
Then follows a detailed explanation of the new method of Lie series illustrated by 
the restricted problem of three bodies. For this a number of comparing calcula- 
tions are made using the power series method and the Runge-Kutta- Fehlberg 
method. This specific method of Lie series compared with W. Grbbner's method 
distinguishes itself by being a numerical method of any high truncation error. A 
simple and efficient automatic step-size control can also be realized. 


1 . The decisive suggestion to use the method of power series as a basis for 
the new Lie series method was made by Dr. Erwin Fehlberg of the Compu- 
tation Laboratory, George C. Marshall Space Flight Center (see also 
Fehlberg and Filippi [ i ] ) . The work was financially supported by means 
of the NASA Contract NAS 8-11209 to the General Electric Company. The 
numerical computations were performed on an IBM 7094, Model II com- 
puter at the George C . Marshall Space Flight Center. 

2. Dozent at the Technische Hochschule Aachen, Aachen, 

West Germany, 



I. INTRODUCTION 


During the last few years the method of Lie series developed by W. Grobner 
[2] was often applied to the numerical solution of problems of initial value in ordi- 
nary differential equations [3, 4, 5]. In this form the method of Lie series could 
in no way compete with other methods for the numerical solution of ordinary dif- 
ferential equations. This was mainly because of the difficulties in finding at first 
a satisfactory analytical initial approximation for the problem considered, and 
secondly because the perturbation integrals to be computed in using quadrature 
formulas required a very considerable computational effort. Furthermore the 
method of Lie series as originally developed had a relatively low order of the 
truncation error, i. e. , the highest order of the truncation error was 0 (h T ) [4] . 

And it was impossible to enlarge the order of the truncation error by adding fur- 
ther terms of Lie series, because the appropriate operator expressions D v 
for v = 4,5,. . . became much too long, and in comparison with the solution the 
effort was by far too much. The method of Lie series described in this paper 
(see also [1] ) eliminates all these difficulties. Now a solution of problems of 
initial value of any high truncation error can be realized by means of a numerical 
method; also included is an efficient automatic step-size control. Only by this 
treatment did the new method turn out to be a very effective numerical method. 

In the first part of the paper the most important theorems of the method 
of Lie series will partly be stated and partly proved. Then follows the new method 
as applied to the restricted problem of three bodies. And finally a number of 
comparing calculations are made with reference to the method of power series 
and the one of Runge-Kutta-Fehlberg. 

II. GROBNER'S METHOD OF LIE SERIES 
A. General 


In 1960 W. Grobner developed a new method for the numerical solution of 
problems of initial value in ordinary differential equations, by using the Lie series. 
During the last few years Grobner' s method sometimes was applied to the numer- 
ical solution of multibody problems. 

In the first part of our paper the most important theoretical basic theorems 
of the method of Lie series will partly be explained and partly be proved. 
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We introduce a linear differential operator 


D = fl(z)_ 9^ +fz(z) 


dz? 


+ f (z) 
n 


JL 

3z 

n 


(i) 


with the holomorphic functions f^( z) of the complex variables z 1( z 2 , . . . , z^. 

We call a function f (z) of these variables holomorphic at one point (e. g. , at 
z t = z 2 = . . . = = 0) if it can be developed as a regular absolutely convergent 

power series. 

If ail f^ (z) and f(z) in the neighborhood of this point are holomorphic, 
then the operator, equation ( 1) , can be applied to f ( z) 


Df(z) = fi(z) 


MLH 

3z t 


+ f2(z) 


9f ( .z). 

dz 2 


+ . 


+ '»<*> 


9z 

n 


(2) 


and we obtain, according to well known theorems of function theory, a function 
that is again holomorphic at the same point. By repeated use of the operator 
(1) tof(z) 

D k f(z) = D [ D k_1 f( z) ] (k =0, i, 2, . . . ) (3) 

we always get a function being holomorphic at the same point. 

Using operator (1) we define a Lie series by 

f -£-D k f( z ) = e'V) (4) 

k =o 

where e*^f ( z) is a symbolic way of writing series ( 4) . 

In equation (4) , t stands for a new complex variable, independent of the 
variables z lt . . . , z^. Every term in equation (4) is a well-defined holo- 
morphic function of these n + 1 complex variables. 

If it can be shown that a number T > 0 exists so that series (4) is abso- 
lutely convergent for 1 1 1 < T, then we can say that series (4) is a holomorphic 

function of the variables z t Zo, . . . , z , t. 

n 
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B. Proof of Convergence 


Using Cauchy's method of majorants [2] we are going to prove now the 
convergence of the Lie series, i. e. , we shall prove the following theorem: 

Theorem 1 


If G is the common holomorphic region of f (z) and of all f^fz) > then a * 
every point of G there is a number T > 0 so that the Lie series 

e tD f(z) = Z "F ° k f< z ) ( 5 ) 

k=o 

assuming 


D = f^z) 


8zj 


+ f 2 (z) 


0z 2 


+ f (Z) — — 

n 8z 


n 


is absolutely convergent for all |t| < T. 

The proof of this theorem, using the Cauchy majorant criterion, is given 
by the following three steps: 

(a) Determination of a suitable convergent majorant for f(z) and for 
operator D, sothatf(z)< g(y) and |D|< A 

(b) Proof that |D^f(z)| <A^G(y) follows from (a) 

( c) Determination of a number T > 0 for which there is absolute conver- 
gence for all 1 1 | < T in series ( 5) 

Step a 

If P is a point within G, then the function f(z) , assumed to be holomorphic 
in G, can be developed into a power series convergent at P. If by an appropriate 
coordinate transformation, P is moved to the origin z t = z 2 = . . . = z^ =0, 
then the power series expansion is 


4 



Q . „ . oo 

f(z) = ) a. , z} 1 . . . z> (6) 

L ii. . . i n 

ii. . . i n 

1 n 

with |zj < p ^ > 0. Power series (6) is absolutely convergent in a certain region 

of the origin. The bounds p are chosen so that convergence of series ( 6) for 

| Zj | < p^ is certain. Generally, it is very difficult to state the whole domain 

of convergence of equation ( 6) . But it is sufficient to find only one domain of 
convergence of equation ( 6) . 

The convergence of 



a 

h 




(V) 


follows from equation ( 6) . Since the terms in a convergent series are bounded 


a. . p. 1 . . . p n < M 

ii. . . i i n 

1 n 


or 


a 

ii. . 




M 

. p 

n 



we get as the majorant for f ( z) 


( 8 ) 


If (z) I ^ g(y i ) = 



b. 

4- 




with |z. | < y. < p.. It follows from equation (7) that 


(9) 




^O. » . °o 



M 

• -< n 


yi 


■/ ln = m y 

n L 


o. ♦ 


JLl 


Jji 


in 


*1 1 „ 

ii. . . i Pi • • • P n n 
1 n n 



(10) 
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because y. / p^ < 1. 

Since an absolutely convergent series within its domain of convergence 
can be rearranged arbitrarily, we arrange equation (10) in the following way 




J b. 

L . : : — b* 

it + i 2 +. . . + i n - 1 



( 11 ) 


the summation i t + i 2 + . . . + i n = i being extended over all possible combinations 

of indices. If we then replace y. and/ or p. with y and/or p so that y = max (y.) 

J J 1 

and/or p = min ( pj) , we get for the function g(y.) with n real variables y. 

function G(y) with one real variable y. So we have 


| f (z) | < G (y) = ) bvy 1 . 


( 12 ) 


l = o 


It follows from 


} |a | (i = 0, 1, 2, . . . ) 

1 ll • • • 1 

ii + ..« + i=i n 

1 n 

and 

Iz.l ^ y.<y < p ^ p. 

that 




< 1 and 


1- 


y i 


1 



(12a) 
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We get 


|f (z) |< g(y.) =£ G(y) = 


(■-f-y 


Since an absolutely convergent series can be differentiated term by term 
as often as desired within the domain of convergence, and the majorizing operation 
can be transferred to the derivatives, we get 


a f(z) 

ag(y k ) 

9G (y ) 

9z k" 

^k 

ay 

8 n f{z) 

< 

a n g(y k ) 

’ a n G(y) 

•v ' 

< 



In the same way as for f(z) we find a majorant for the functions f ( z) 

K 


iyz) i ^ 


(■ -?r 


We assume that the p. are equal for all estimates; this can, of course, always be 
achieved. 


We therefore get as majorizing operator 
. N d 


I D| < A = 


(-rf 


. ( N = k max {N^} ) 


It follows from 
I f(z)| ^ G(y) and | D | < A 


|D k f(z) | < A k G(y) . 



power series 


Since the functions f^( z) are holomorphic, they can be expanded into a 


V z >= I 


O. . . oo 


(k) 


C.' . Zi Jl . . . z Jn 

r Ji •••]„ n 

Jl. n 


(18) 


We also develop the majorizing operator in a power series 




dy 


j = o 


(19) 


With f( y) > IfjJ 2 ) l» follows that 


Jl 


(k) 

c i 

ii* • • J 

b +. . . + j = j n 

J n J 


• (j = 0, 1, 2, . . . ) 


(20) 


Assuming this, it now follows that 


Df ( z)= J c< k > . z^ 1 . . . z d ^ a. . z/ 1 . . . z 1 

z Jl- - - J n n dz k L ii- - - i n n 


Jn d_ 




z l 


h + Ji 


, % + jk - 1 
k 



n 


(21) 


and 


AG (y) = ) i y.b.y 1 + j 1 . (22) 

With equations (12a), (20), and i < i (i t + i 2 + . . . + i = i), these 
relations hold: n 

1. 

2 . 




ZZ I 


I 


k = 1 i t +. . . +i = i 
n 


Ji + - • . + J n = j 


, 0 < k > 

Jl* • 


n 


a. . l 

ii« # .1 
1 n 


1^ y b. 
J i 
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1. e. , AG(y) is a true majorant for Df(z) because the coefficients of the series 

2. are always smaller than the corresponding coefficients of series (22). Thus 


I Df ( z) | < A G (y) . (23) 

But since Df(z) and/or A G(y) , as was shown above, fulfill the same conditions 
as f ( z) and/or G(y) , we get 


D[Df(z)] <A [A G(y) ] 


D [D k_1 f(z)] < A[A k-1 G(y)] 


(24) 


Since the sum of the majorants is majorant for the sum of the corresponding 
mino rants, we get 

I IZ — D f (z) | < LZ^A G(y) . (25) 

k=o k=o 


Step c, 

To determine the convergence radius of the majorant Lie series (25) , 
we first calculate from 
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The proof of equations (26) is obtained by mathematical induction. The validity 
of equations (26) for k= 1 has already been demonstrated above. 

From equations (25) and (26) it now follows that 


OO , & I 


k -o 




Series (27), and also the Lie series (25), converge for 


-N (n + 1 ) t 
n+ l 




< i , 


i.e.,for 


It! < 


(n + i)N 


<rt) 


n +i 


= t 


(28) 


Thus N depends on operator D, while p depends on f( z) as well as on D, and y 
can be freely selected between 0 and p . At the expansion point P(y =0) the 
series is absolutely convergent for 


1 1 1 < T* = 


(n + 1) N 


By comparison with the binomial series 


(29) 


} <~k> xk= < 1 + x >~ k ( 30 > 

k =o 


we obtain an evaluation for the sum of the Lie series by taking as upper bound 
the summation value for the majorant 


I 


k =o 


t k k 

l^ Df(z) 




-N (n + 1) 




n 

n + i 


( 31 ) 
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E stimate of Remainder 


To estimate the error of the Lie series 

) ” ~~D k f (z) = e tD f(z) (32) 

k = o 

when it is truncated after the m-th term, we get from equation (27) 



since 


(33) 


ilL 

T 


< 1 


NOTE: If the f^(z) are polynomials, the estimate of the remainder (33) 
can be further improved. 


ii 



C. The Lie Series for the Solution 
of Ordinary Differential Equations 


The following important theorem [10] can be derived from the analyses 

in n-B: 


Theorem 2 


If G is a finite, closed region in z-space where all f^ (z) and f(z) are 

holomorphic, then there exists a positive number T so that Lie series (4) in 
II-A is absolutely and uniformly convergent in the entire region G for |t| < T 
and is a holomorphic function of the n + 1 complex variables z 1; z^ . . . , z^, t. 

Lie series ( 4) can be differentiated term by term as often as desired 
within G and |t | ^ T with reference to the variables z 1# z 2 , . . . , z n , t. Thus 


and 



(34) 


(35) 


NOTE: In the special case f (z) = Z. we get the Lie series 

Z.=e tD z. . (i = i, 2, . . . , n) (36) 

In a certain region of the variables z t , Z 2 , . . . z in which all f (z) of 

n k 

operator (1) in II-A are holomorphic, Lie series (36) represents, according to 
Theorem 2, for |t I < T, holomorphic functions of z lt z 2 , . . . , z^, t. For 
t = 0 the functions Z,, in equation (36) take the initial values z^ 

(Z.) = z . (i =1, 2, .... n) (37) 

1 t = o 1 

The following important properties can be demonstrated [2] for the Lie 
series defined by (2) in II-A: 
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( a) For the sums and products of Lie series formed with the same 
operator D: 

et ° [cigi(z) + c 2 g 2 (z) + . . . + c m & m ( z )l 

= c 1 e tD g 1 ( z) + c 2 e tD g 2 (z) + . . . + c m etD g m ( z ) ( 38 ) 

where the c are constant, and 
1 

e tD [gi (z)g 2 (z) ♦ • • S m ( z ) 1 = fe tD gi(z) ] [ e tD g 2 (z) ] . . . [e tD g m (z) ] . 

(39) 

(b) Theorem of permutability: 

If F(z) is an arbitrary holomorphic function near z lt z 2 , . . . , z^ 

and its power series expansion also converges at the point {z ls z 2 z } , 

then the functional symbol F and the symbol e^ for the Lie series 


F( Z) = £ 


k = o 


t k 

TT D f(z) 


or 


(40) 


_ . tD . tD_, . 

F (e z) = e F(z) 

can alternate. 

Also, according to equation ( 34) , the special Lie series in (36) can be 
differentiated term by term beyond t in the region 1 1 1 < T. Thus we have 


9Z. 


Z.= 

l 


at 


1 


k = 1 

°o k 


kt k_1 n k y-^-t k k+1 

D z. = / ~ D z. 

l *— 


k(k-l) ! 


k = o 


k! 


“ _t-k. 


= EZ ir D <D2 i ) = ntr Df i (z)= V z > 


k = o 


k = o 


(i = 1, 2, . . . , n ) 


(40a) 
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The functions 1 ( z) are holomorphic in the neighborhood of the initial point zj, z 2 
. . . , z n and the independent variable t is not explicit. 

According to ( 40a) the functions Z . in ( 36) satisfy an autonomous system 
of first-order differential equations. From this we immediately obtain the theorem: 

Theorem 3 

K we have an autonomous system of first-order differential equations of 
the form 

Z = f.(Zi, Z 2 Z ) (i = 1, 2, . . . , n) (41) 

li n 

with the initial conditions 

<Z.) t=o = z. (i = 1, 2, . . . , n) (42) 

and in whose neighborhood the functions f. (z) are holomorphic, then the Lie 
series 



where 


D=f * (z)- !r Mz) "4 + + f n <z) 


n 


represent the desired solutions of this system of differential equations. 


Thus the z. should be considered primarily as variables to which operator 

v 

D is applied. After all the differentiations, prescribed by the operator D z^, 
have been performed, the z. can be replaced by the initial values. 

Supplement 


A nonautonomous system of first-order differential equations 

Z -f.(Z lt Z 2 Z , t) (i = 1, 2, . . . , n) (44) 

ii n 

can always, by the substitution 


t = 


Z 


n+i 


(45) 
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be converted into the autonomous system 


Z.= f (Z t , Z 2 , . . . , Z n , Z n+1 ) . (i=l, 2, . . . , n, n+i) 
It follows from equation ( 45) that 


dZ 

n+i _ dt _ ^ = £ 

dt dt “n+i 

and ( Z ) = t 0 . 

n+i t=t 0 


( 46 ) 


In general, Lie series ( 43) converge very slowly and are therefore not 
usable in this form in a numerical computation. This disadvantage can be elim- 
inated by rearranging (43). To achieve this, operator D is split into 

D = Dt + P 2 (47) 


where 


n 


D i - L_<P j ( z l» z 2» « ♦ •, z p ) 9Z 
i=l 


and 


n 


Dfe= 7 [f.tZi, Z 2 , . . . , Z n ) - cp.iZi, Z, 2 ,.. . , Z n )] — . (48) 
i = 1 1 

Thus near the initial point | L - cp A < I <p . I and the functions belonging to Dj 


oO 

z. (t) = } ^rr o) [D^Zil . . (i = i, 2, . . . , n) (49) 

ia k = o Z (0) 

should be regular and expressed by known functions in closed form for all finite 
values of t. 

The subscript Z ^ ^following the brackets in equation ( 49) indicates that 

the variables Z u Z 2 Z are to be replaced by the initial values (42) only 

n k 

after all the differentiations prescribed by D have been performed. 
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By splitting (47) the desired solutions Z^(t) can be written in the following 
form ( see Grobner [2] for details) 


V‘ )= Z ia (t) + l f 


t .. .m 

(t - t) 


m_ 


m = O tn 


, [D 2 D Z.l 
ml 1 1 1 Z 


dT. 


ia( T ) 

(i = 1, 2, . . . , n) (50) 


NOTE: If the term 9/9t is added to operator D, the splitting of operator 
D becomes considerably more flexible with respect to possible 
choices of the functions <p\. By splitting D into 


Di 



<P i ( Z l> Z 2» 


n 


t) 


9 9 

9 Z . + 9t 
1 


(51) 


and 

1 ^ - /_ Z 2 > • • • > z n ) - cp.(Z u z 2 , 

i = 1 


* * Z n’ t ^dZ 

1 

(52) 


we can choose any desired functions of t — polynomials, trigono- 
metric sums, m-terms of the power series expansions, etc. — as 
functions <p 


D. Generalized Lie Series 


To extend the scope of Grobner's Lie series, Filatov [6] introduced 
"generalized Lie series": Under the assumptions that we have <^(1* z i^ 2 > • • • » z n ) 

(m = 0, 1, 2, ..., n) holomorphic functions of the complex variables t, z l5 Z 2 , 

. . . , z in a region G of the (n + 1) -dimensional space R ( t, z t , z 2 , . . . , z ) 
n n +1 n 

containing the origin of the coordinate system. Furthermore, the function 

f(t, Zi z 2 , . .. , z ) may also be a holomorphic function in the same region G. 
n 

We now introduce operator W 
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9 Q 

W = <po(^> z l z 2» . • • , Z^) ^ ■*" 2 k^’ z l* Z 2> • • • ’ z n ) g z ( ^3) 

k = 1 k 

for which these relations apply when operator W is repeatedly applied to function f: 
W° = f and W n f = W ( W n_1 f) . ( 54) 

We can now, just as in equation (4) , define the "generalized Lie series" 

oo k 

'S 4 „ 7 k, ,. , tW 

2 -£T" W f (t, z t , z n ) = e f. (55) 

k = o 

From equation (55) we again get Grobner's Lie series, for the special 
case (p o=0 and for the functions <p (k = 1, 2 n) and f which are indepen- 

dent of t. 

Absolute and uniform convergence in a finite and closed region G of R - 
space can also be demonstrated for these "generalized Lie series" [6] . 


III. APPLICATION TO THE RESTRICTED 
PROBLEM OF THREE BODIES 

A. General 


Since the Lie series for a given initial value problem generally converge 
very slowly, it is usually necessary to split operator D into two operators 
D = Dj + E >2 and use formula ( 50) for the calculation. But in doing this it is 
necessary that the functions belonging to Dj [see equation (49) ] represent a 
good approximation of the desired solution, and have to be regular for all finite 
values of t, and can be expressed by known functions in closed form. Until now, 
a low-order polynomial, second or third-order, in (t - t 0 ) with the coefficients 
D z. has been used at each initial point. 

The operators D^z. for v = 0, 1, 2, ... had to be computed in this way, 

but generally this is possible only up to v = 2, 3, or at most 4, because the oper- 
ator terms soon become too long. A good initial approximation was rarely 


17 



achieved and the highest order of the error term for the Lie series method was 
very low [4] . 

Since the Lie series for a problem in differential equations, for reasons 
of the uniqueness of the solution, is simply the Taylor series for the solution 
function, the power series method can obviously be used to get a good initial 
approximation for the method of Lie series [2] . Thus it becomes possible to 
expand the method of Lie series into a numerical method of any desired order 
of accuracy. By adding perturbation integrals to formula (50) , the order of 
accuracy of the method of Lie series can be raised by a full h-power over the 
power series method by every further perturbation integral. The perturbation 
integrals are therefore analogous to the k-values of the Runge-Kutta-Fehlberg 
method. 


B. The Restricted Problem of Three Bodies 


For the numerical solution of the restricted problem of three bodies in a 
rotating coordinate system the problem of initial value [7, 8, 4] 




& 2 


x + 2S 2 - M**' 

X + /i 

[(x + m) 2+ y 2 ] 3/12 P 

y - 2d-\ - li ' 

y , .. 

[(x+m ) 2 +y 2 ] 3/2 M 


x - /u 


;Tyj 37T = ^3 (56) 

X 


with the initial conditions 


x (0) = x 0 , 

x ( 0 ) 

= x 0 

y (o) = y„ 

y (0) 

= yo 


(56a) 


will be integrated step by step. Additionally, n = 1/82. 45 and 
With the abbreviations 


^ [( X +H) 2 + y2] 3 /2 = Nl 


= 1 - M. 


(57) 


“ I(x-m ') 2 +y 2 ] 3 / 2 = N 2 
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equations (56) become 


x = 

y = ^2 

$1= = X + 2^ 2 - ^ (x +M) - ^ (x - M ') (58) 

^ ^ = y-«i- ^ . 

Operator D [ see equations (43) or (51) ] for equations (58) is 


_a_ 

dx 


_a_ 

oy 


D at + ** a x + ’ is dv + a^< + * 


<V 2 


(59) 


We split operator D in equation (59) into 


^ _ „ a „ a a d a 

Dl “ 1?1 ' ax + ^ 2 »" + (pz ~^~ + ( Pi~ T ~ + 


oy 


a^! 


9t?o at 


and 


° 2 " + <^4-<?4) 


(60) 


where 0 = 0!+ D^, and <p 3 and are arbitrary functions that will later be 
chosen in a more convenient form. Using equation (50) , the solutions of equations 
( 58) generally take the form 


oo t k 

f(t,z)= g(t,z)+ ) J (t k ~ S) [D^z] .ds 

/N * ' 


k = o o 


g) 


(61) 


where g (t, z ) = e^z must be a known approximate solution of equations (58) . 
We now choose cp 3 and as follows: 

(P3 = ltfx] 0 +...+ [D m+2 x ] 0 

III • 

r _«> , ( t -t n ) m r „m +2 , 

</ , 4 =: [D 2 y]o +...+ I D y ) o 


(62) 
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V V 

where [D x ] 0 and [D y ] 0 represent the v-th total derivatives of x and y with 
respect to t, taken at the point t= t 0 , x= x 0 , x = x 0 , y = y 0 , and y = y 0 . 


With <p 3 and <p 4 chosen according to (62) , the approximate solutions g(t, z) 
of (62) are 


x a = [e tDl x ] 0 
y a = [e tDl y ] 0 
x = [ e tDl ^!] 0 

a 

y a = [e tDl ^23o 


(63) 


where [ ] 0 means that the expression in brackets, after having performed all 
the operations prescribed by D 1; is computed at the point t = t 0 , x = x 0 , x = x 0 , 
y =y 0 , and y = y 0 . If, for example, we calculate x from equations (63) we 
get, using 


D^X = X 


DjX = d - 1 

Di 2 x = D 4 ( D lX ) = =<p 3 = [tfx ] 0 + . . . + ^-iD m+ 2 x ] 0 

D t 3 x = D t ( D t 2 x) = D^s { [tfx] o + . . . + [D m+ 2 x] 0 } 


= [D^xlo + 2 gV ^ [D 4 x ] 0 + . ♦ . +m 


\ , u t 

m! 




Dj 4 x = Dj ( Dj 3 x) = -J-r { [D ? x] 0 + -^T^tl^xlo + . . .} 


2 ! 


m -2 


= ID‘x] 0+ ^Vx ] 0 + . . . + “USlUltkJ lD m+ 2 x] 


_ m r m +2 , 

D l x = [D x ] 0 

Di m+;l x = 0 for j = 1, 2, 3, ... 
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under consideration that in x = [e tDl x ] 0 , D^x has to be taken at the point t = t 0 , 

3. 

the following equation for the approximate solution x : 

ct 


X 

a 



[D^xlo 



[ D^x ] 0 


In the same way, we can show that 


a 


= ) ~ [Dj^j] o = > m 2 [ d ^!] 0 

v\ L v\ 

ii = n 7» — A 


V = o 

oo .... V 


v=0 

m+2 v 


y a = T [Di”yl 0 - r~ (t ,: M I pl, y] 


V = o 


V = o 


y = o 

°° /JL . m+2 , , , v V 

k - . 

V = o ' 


(64) 


NOTE: The solution of the system of differential equations 


dZ. 

-^ _ =t?.(z 1 , Z 2 , . . . , z ) (65) 

with 

Z.(t 0 ) = z. (0) (i= 1, 2, . . . , n) 

can, using the formula of Lie series, be written in the following 
closed form [see equation (43)] : 


: i<*> = I 


oo V 

ilztfl) V 1 

v\ 1 i J z(0) 


v= o 


Operator D is here defined by 


D= ) ^k (Zl ’ Z2, * * * ’ Z n ) 


k=l 


9z, 


(66) 


(67) 
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Given equation (65) , i. e. , given 


dZ. 
1 

dt 



9 


it follows from equation ( 67) that 


Dz. 

l 



dZ. 

l 

dt 


And from equation (68) we get 


T?z. = D ( Dz.) = IW . 
i i l 




dj. 

*k~s r 

k 


Sz. 


( 68 ) 


(69) 


so that Lie series ( 66) , as is also shown by the uniqueness of the solution, is 
identical with the Taylor series for z^(t) . 

Therefore it is possible to get the desired accuracy for the initial approx- 
imations (64) by any choice of m using the method of power series [7] without 

p 

explicitly calculating the operators D z. . The explicit calculation of the oper- 

V 1 

ators D z^ for v = 4, 5, . . . leads to extremely long expressions being of no 

value to the practical solution [4] . It was mainly for this reason that the method 
of Lie series could not compete with other methods for the solution of ordinary 
differential equations. By use of the method of power series, however, these 
difficulties could be eliminated. Furthermore it is now possible to expand the 
method of Lie series by making it a numerical method with any high truncation 
error possible. This was the decisive step in the new method of Lie series first 
developed by Fehlberg and Filippi [ 1 ] . 

Parallel to the k-values in the Runge-Kutta-Fehlberg method [1, 8, 9, 10], 
we add perturbation integrals [see equation (61) ] to the initial approximations to 
raise the order of the truncation error. To solve these perturbation integrals, 
k 

the operators t^D z. for k = 0, 1,2,... must first be computed. Each pertur- 
bation integral added to the initial approximation (64) enlarges the order of the 
truncation error 0(h) by one full power. 
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The following equations are the results of our calculations for k = 1, 2, 3, 

and 4. 


(a) k= 0: 

(b) k= 1: 

(c) k= 2: 

(d) k= 3: 


D>D°x = 0 
E^D°y= 0 

DjDx = t> 3 - <ps 
DzDy =^4-^4 

D 2 D ? x= 2(^4 - (pi ) 

E^tfy= -2(t 9- 3 - <p 3 ) 

D 2 D 3 x= U 3 -<p 3 )(I) + (^4 - ^»4>< 11 ) 
Dzrfy= u 3 -n)(iD + ^ 4 -<P4)(in) 


with 


(I) = - 3 - M ' 


y 2 - 2 (x+M) 2 
[ (x + a 0 2+ y 2 J 5 / 2 


y 2 - 2 (x - /Q 2 

f(x-M') 2 + y 2 ] 5/2 


(ID 


m' 


3y ( x + m ) 

[(x + M ) 2 +V ] 5/2 


+ 


3y(x-/Q 

[(xVr + y 2 ] 5/2 


(HI) 


= - 3 - m' 


(x + M) 2 - 2y 2 

[(x + M) 2 + y 2 ] 5 / 2 


- M 


(x - /-Q 2 - 2.V 2 
( (x-/x ) 4 + y 2 ]V 2 


(e) 


k= 4: D2D 4 x= 2(^ -(p 3 ) [^i(IV) + i? 2 (V)] 

+ 2^ -<p 4 ) [o?i(V) + t? 2 (VI) + (I) + (HI) + 4] 
2(^4- ^ 4 ) [ ^i(VI) + ^ 2 (vn)] 

+ 2(^3 -<^ 3 )[^i(V) + tMVI) - (I) - (in) - 4] 

with 


<IV) ■ T **V$*+W [3y! - 2(xt><)21 + [(x-V'tC^v^ 3 ^- 2 ^-^) 2 ] 
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(V) = 

3yM' 

[y 2- 4 (x + M) 2 ] + 

3 y^ 

[ ( x + fi) 2 + y 2 ] 7 / 2 

[ (x - ju') 2 + y 2 ] 7 / 2 

(VI) = 

3 ( x + 

[(x + M ) 2 - 4y 2 ] + 

3 (x - m')M 

[ (x + n ) 2 + y 2 ] 7 / 2 

[(x - n'f + y 2 ] 7 / 2 

(VII) = 

3yn' 

■ [3(x+ M) 2 -2y 2 ] + 

3ny 

[(x + nf + y 2 ] 7 / 2 

[(x-/i') 2 + y 2 l 7 / 2 


[y*-4(x-M') 2 ] 


[(x-m') 2 - 4y 2 ] 


[3(x-M') 2 -2y 2 ] 


We shall now combine the solution formulas for equations (58) , consider- 
ing also equation (61) with the perturbation integrals for k = 0, 1, 2, 3, and 4: 

Solution Formulas 

t t 

x(t)= x a (t) + f (t - t) (i> 3 -<p 3 ) a d r + / (t - t) 2 (^ 4 -(p 4 ) a dT 

to t 0 

1 } 

+ J J T) 8 [(* 3 -«P3)d) + (*4 -^4)OD] a dT 

t 0 

t t 

+ ^ / (t- T ) 4 {(<h-<p 3 )Ui(IV)+ * 2 (V)] + (i>4-^4)[^l(V)+d2(VD 


+ (I) + (in) + 4]} dr 


a 


t t 

y(t)= y (t)+ / (t - t)(i>4 -(p^ dr - / (t - t ) 2 (^ 3 -<^ 3 ) a d T 
t 0 t 0 

+ T J (t-r) 3 [(^ 3 -<^ 3 )(n)+ (i> 4 -<p 4 )(IID] a dT 
to 

! t 

+ U / <t-T) 4 (U4-<?4)Ui(VD+ tfi(vn)] 


+ (* 3 -<Ps) Ui(V)+ 1? 2 (VI) - (I) - (III) - 4]) a dr 
t t 

x(t)= x & (t) + f(f 3 -<p 3 ) a dr + 2 f (t -T )(^4"^4) a dT 
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+ J- f (t-T ) 2 [U 3 -<p 3 )d) + U 4 ~<Pd (ID ] dr 
to 

t t 

+ ~rf (t " T > 3 ^3 - ^3) Wi(iv) + ^ 2 ( V )] 

to 


+ (t> 4 - <P 4 > Ui(V) + ^ 2 (VI) + (D + (III) +4]} dr 

d. 


t t 

y(t)= y a (t) + f -<Pi) d t- 2 f (t - t) ( i >3 - cPs) a dr 

t 0 t 0 


1 t 

+ 7 - / (t - r ) 2 [(^ 3 -^> 3 ) (II) + (* 4 -? 4 ) (HD! dr 


1 1 

+ — / (t *r ) 3 {(t ?4 - 1 ^) 4 ) [i?i( VI) + ^ 2 (VH) ] 


+ ( 1 ? 3 - <p 3 ) [ ^ !( V) + * 2 (VI) - (I) - (III) - 4]} dr 

d 


Step-by-Step Arithmetical Performance, 


By use of the solution formulas in the second part of this section, we can 
immediately perform the step-by-step computation of x(t) , x(t) , y(t) , and y(t) 
if we compute the perturbation integrals that depend on m. The m is arbitrarily 
chosen for computation of the approximate solutions x ( t) , x ( t) , y ( t) , and 

y^(t) via a quadrature formula with an error of the order 0( h ) , 0(h ) , 

0(h m +3 ) , or 0(h m+ 4 ) . But the effort is by far too great [4] . 


We shall now suggest a more satisfactory method for calculating the 
perturbation integrals. 
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Calculation of the First Three Perturbation Integrals 


We start as follows: 


Us - (« - t») m+1 A m+1 +(t -**> m+2A m+ 2 + (t - *• > m+3A m + 3 < 70 > 

<’’* - ^>a = (t ' W m+ ‘ B m + l + (t - t » >m+2B m +2 +(t -‘» )m+3B m + 3 . < 71 > 


If we successively introduce t=t 0 + 2A t, t =t 0 + A t, and t = t 0 - A t into ( 70) , we 
get 


A . +2A tA +4(At) 2 A = 

m+1 m+2 m+3 


A , + At A + (A t) 2 A 

m+1 m+2 m+3 


(t? 3 - (fi 3 ) 


t=t 0 + 2At 
a 


(2A t) 


m+1 


= a 


(^3 ~<P 3 ), 


t=t 0 + At 


(At) 


m 


+1 


-= fi 


(72) 


(73) 


t =t 0 - At 

(^3 - ^3) 

A - AtA + (A t) 2 A _ = — — — 

m+1 m+2 ’ m+3 . , , m +1 . A . . m+ 1 

(-1) (At) 


'= y 


(74) 


We get three corresponding relations if we successively introduce t = to + 2At, 
t = t 0 + A t, and t = t 0 - A t into equation ( 71) . 


If we add equations (73) and (74) , we get 


(75) 


If we subtract equation ( 74) from equation ( 73) , we get 


2AtA m + 2 “O' V 


(76) 


or 


A = ft ~ Y 
m+2 2A t 


(76 a) 


Using equation (76a) we get from equation (72) 
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or 


+ 4(At) 2 A 


= o; - /3 + y 


m+1 ' ' m+3 

If we subtract equation (75 from (77) , we get 


(77) 


3 (At) 2 A 


m+3 


!£ + jl 

2 2 


m+3 


_ 1 / 3/3 y\ 

3 (At) 2 “ 2 2 J 


(78) 


(78a) 


It follows from equation ( 75) , using equation (78), that 


a = -2L 
m+ 13^3 


(79) 


We then get these expressions for A^, A m+ 2 ’ 


and A 


m+3' 


m + l (it) m + l ' 6 


1 t =t 0 + 3A t t = to + At 

~(^3 “^3) q -6(i> 3 - (pz) 

III <x a 


m 


+ (-l) ^(^-^g) 


2 

t =to - At 


(80) 


m+2 (A t) m+2 


^-[(^3-<P3)* t0+At +(-l) m (^ -<p 3 ) 


t =t 0 -At 
a 


(81) 


m+3 


m+3 


+ (-l) 


(At) 
m+ 1 


± 

6 


1 . t = t 0 + 2 At t =t 0 + At 

~ ( $ 3 “ <P3)„ -3(i? 3 - <p 3 ) 

in d a 

z 


(^ ~<P3) 


t =tfi - At 


(82) 


Analogous to this we get three expressions for B ® m +2’ anc * B m+3' 

We simply say B m+1 instead of A ^ etc., in equations (80) , (81), and (82), 

and we change every (j? 3 - (pz ) to (i ?4 - (Pi) • 

a a 

If we now introduce the expressions ( 80) , ( 81) , ( 82) , and the corre- 
sponding terms B „ , B and B , _ to the first three perturbation integrals 
m+1 m+2 m+3 
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we get, by means of a term-by-term integration, the following definitive solution 
formulas: 


x(t) = x (t) + 

a 


1X1+3 A ....m +4 A ... . 111+5 

A .(At) A C (A t) 

+ m+ 2 t } m+ 3 

(m+2)(m+3) (m+3)(m+4) |(m+4)(m+5) 


A ..(At) 
m +1 


2B m+1 


m+4 


m+ 5 


(m+2 )(m+3) (m+4) 


(m + 2) (m+3) (m + 4)(m + 5) 
m+ 2 


2B „ (At) 
m + 2 v ' 

(m + 3) (m + 4)(m + 5) 
m + 5 


(83) 


x(t) = x (t) + 

a 


2B m+1 < At > 


A m +i (At) 

m+2 

m+3 


A m+ 2 (At > 


m+ 3 


A m+ 3< At > 


m+ 4 


(m + 2) (m + 3 ) 


m+3 

2B m + 2 < At > m * 4 

(m + 3 ) (m + 4) 

m+ 4 


m+4 


(84) 


l(I)A m +1 + <n > B m + iy At > 


(m + 2) (m + 3) (m +4) 

B m+l (At)m+3 
y(t) y a (t) + ( m + 2 ) (m + 3) 


, B m+2 (At > m+4 1 i B m + 3 (At > 

+ — r— — +i 


m + 5 


2A m + l (At) 


m+4 


( m + 3) ( m + 4) i ( m + 4) ( m + 5) 

2A m+2 < A ‘> m+5 


( m + 2) ( m + 3) ( m + 4) 


» n > A m ,l +(1I1 > B m+ lla< At > 
( m+ 2) (m + 3) (m + 4) (m + 5) 

m + 2 


(m + 3) (m + 4) (m + 5) 
m + 5 


(85) 


y(t)= y (t)+ 

cl 


Vti< A y_ 

m+2 
m + 3 


B .(At) 
m+2 


m+3 


B q(At) 
m+3 


m +4 


2A m+l (At > 
(m + 2) (m + 3) 


2A (At) 
m + 2 v ' 


m+3 
m + 4 


m+4 


(m + 3) ( m + 4) 


( 86 ) 
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(m + 2)(m + 3)(m + 4) 


m+ 4 


where (I) , (II) , ( III) are abbreviations used in the expression for k = 3 in 

Si Si Si 

III-B, d (page 23) . Index a means that in (I), (II), and (HI) x, x,y, y always 

has to be replaced by x , x , y , y [see (64) ]. The expressions A 

a a a a m 1 


!!•••» 


B are found in expressions (80), (81), (82), and in three corresponding 

m r — ~i 

expressions. The sign J ! indicates that this term is used for automatic 

step-size control. All other terms with (At) m+ 5 are to be neglected when solu- 
tions ( 83) , . . . , (86) are computed step-by-step. 

Computation of the First Four Perturbation Integrals in III-B. 

Using the formulation 


<■>»-«>„= <t - t„) m +1 A m +1 +... +(t-t,) m+4 A m+4 


(87) 


m +1 


( ^ 4 ^4)^ ( t to) ^ • • • "^(t to) 


m +4, 


B 


m +4 


and computing as in the third section of III-B, we get as the definitive solution 
formulas 


....m+3 

x(tl= x (t) + A (At) ' +a — (Al) 

a' m+1 (m+2)(m+3) m+ 2 (m + 3)(m +4) 


+ A 


m+ 5 I m+6 j 

(At) + • ^ (At) l 

m+3 (m +4)(m+ 5) J m+ 4 (m+ 5)(m + 6) | 

1 1 


+ 2B 


(At) 


m +4 


+ 2B 


(At) 


m + 5 


m+1 (m+2)(m+3)(m+4) m+2 (m+3)(m+4)(m+5) 


+ 2B 


(At) 


m+6 


m+3 (m+4)(m+ 5) (m+6) 

. m+ 5 
(At) 


(m + 2)(m + 3)(m + 4)(m +5) [f^^m+1 + ^^^m+l]< 
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(At) 


m +6 


( m + 3) ( m + 4) ( m + 5) ( m + 6) I m + 2 

m + 6 


|(I)A„_ + (IDB m+ 2 j a 


+ 2 


(m + 2) (m + 3) (m + 4) (m + 5) (m + 6) I + 1 [*a^^ + ^a^^]a 


+ B 


m 


+ 1 [x a (V )+ y a (VD + (I) + ( HI) + 4J a [ 

y(t) = y (t) + B 


( 88 ) 


_(At) m+3 (At ) m +4 


m+1 (m+2) (m+3) m+ 2 (m+3)(m+4) 


+ B 


B 


lAt)™^ ' 


m+3 (m+4)(m+5) j m + 4 (m+5)(m+6) j 

i j 


- 2A 


(At) 


m + 4 


- 2A 


(At) 


m + 5 


-2A 


m+1 (m+2)(m+3)(m+4) m+2 (m+3)(m+4)(m+5) 


m + 3 ( m+ 4) ( m+ 5) ( m + 6) 
m + 5 


(At) 


( m + 2) ( m + 3) ( m + 4) ( m + 5) ^ ^m + 1 + ^m + 1 ] a 


(89) 


(At) 


m +6 


(m+ 3) (m+ 4) (m+ 5) (m +6) U) A m+ 1 + ( IH) B m + 1 a 


+ 2 


(m +2) (m+ 3) (m+ 4) (m +5) (m+ 6) ( B m + 1 [*a^^ + ^a^ V ^ ]a 


m + 6 


+ A m+1 fx a (V) +yJVI) - (I) - (III) - 4 


] 


m+2 m+3 /A .,m+4 

x ( t) = x ( t) + A -l At) + A — ^ — + A — 

a m+1 m+2 m+2 m+3 m+3 m+4 
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* 


+ A 


(At) 


m + 5 


m + 4 m + 5 


+ 2B 


(At) 


m + 3 


m +1 (m+ 2) ( m + 3) 


+ 2B 


(At) 


m + 4 


m+2 ( m+ 3) (m+ 4) 

iM m+4 


+ 2B 


(At) 


m+ 5 


m+3 (m+4)(m+5) 


(m+ 2) (m + 3) (m+ 4) [ (I)A m+ 1 + (II) B m + l]a 


(At) 


m + 5 


(m + 3) (m+ 4) (m + 5) ['‘'“m + 2 '“'“m + 2 ] a 


(I)A_ +(II)B_ 


(90) 


+ 2 


(m + 2)(m+3)(m + 4)(m + 5) [ A m + 2 [*a^ ^ + ^a ^ ] 


m+ 5 


a 


+ B m + l[ !i a (V)+!? a (VI) + (1 > +(III)+4 ]a} 


, . m+ 2 

y(t) = y (t) + B — 

yK ’ V ' m +1 m+2 


+ B 


1M 


m+ 3 


m+ 2 m+3 


+ B 


1M 


m + 4 


m+3 m+4 


+ B 


_(M 


m + 5 


m+4 m + 5 


2A 


(At) 


m+3 


m+1 (m+2) (m+3) 


2A 


jAt} 


m+ 4 


- 2A 


(At) 


m + 5 


m+2 (m+3) (m+4) m+3 (m+4)(m+5) 

m+4 


(91) 


(At) 


(m+ 2) (m + 3) (m +4) 
m+ 5 


(II) A . . + ( III) B . . 1 
m+ 1 m+ 1 a 


-(At) 


(m + 3) (m+ 4) (m+ 5) |^ I ^ A m+2 ^ BI ^ B m+2ja 


+ 2 


1M 


m+ 5 


(m + 2) (m + 3) (m + 4) (m + 5) ) m + l|“a ' ]a 


{ B m + lp 


i (VD+V 
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+ A m*l[V V > + *a< VI > - (I) - (III) -%} 


with 


A m +1 = “e" < - O' + 4/3 + - 6 ) 


1 


A 

A 


m+2 12At 
1 


m +3 6 (At) 

1 


m +4 12 (At) 


( - a + 8(3 - 8 y + <5 ) 

2 (oi - (3 - y + 6 ) 
3(0- -2(3 +2y - 6) 


where 


(^3 “ ^3) 


t =t ft + 2A t 


(i ? 3 - <p 3 ) 


t =tn - At 


a = 


P = 


( 2 A t) 

(^3 - cp 3 ) 


m +1 


t =t 0 + At 


(A t) 


m + 1 


r = 


<5 = 


(-At) 


m + 1 


(t ? 3 - <p 3 ) 


t =t n - 2 At 


( -2 A t) 


m + 1 


and also corresponding expressions for + + B m + and + 

where h ? 4 - ip 4 ) appears in the expressions for a , (3 , y, and 6 instead of (^ 3 - <ps) . 

3. ** 

The expressions (I) , (II) , . . . , (VII) appear in III-B, d and e (page 23) for 
k = 3 and k = 4. All other expressions correspond with those in the third section 
of m-B. 


Some Numerical Results 


The following table (see also Fehlberg and Filippi [1] ) illustrates the 
efficiency of the new Lie-series method using three perturbation integrals and 
automatic step-size control as described in the third part of III-B. The restricted 
problem of three bodies [see (56) ] was solved for the initial conditions 

x 0 = 1.200000000000000 

x 0 = 0 

(92) 

y 0 = -1.04935 75098 30320 
yo = 0 


32 


9 


For comparison, the problem was also computed using other modern methods of 
numerical integration: the method of power series and two versions of the Runge- 
Kutta-Fehlberg method [9, 10] . All four methods were applied with a truncation 
error of the order 0(h 13 ) . All computations were performed on an IBM 7094, 
Model II computer in double precision (16 digits) . 


NUMBER OF INTEGRATION STEPS AND COMPUTER RUNNING 
TIME AFTER 12 REVOLUTIONS (ABOUT ONE YEAR ACTUAL TIME) 


Method 

Number of Steps 

Running Time fmin] 




PSE * ) 

5896 

2.49 

Lie * ) 

5191 

2.05 

RKFe 1964* ) 

4740 

1. 83 

RKFe 1965 * ) 

3353 

1.35 


PSE = Power Series Expansion 

Lie = New Lie-Series Method with Three Perturbation Integrals 
RKFe 1964 = Runge-Kutta-Fehlberg Method 1964 
RKFe 1965 = Runge-Kutta-Fehlberg Method 1965 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 
Huntsville, Alabama, August 29, 1966 
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"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof.” 

— National Aeronautics and Space Act of 1958 
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